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Abstract. We examine the effect of accuracy of high-order spectral element methods, 
with or without adaptive mesh refinement (AMR), in the context of a classical 
configuration of magnetic reconnection in two space dimensions, the so-called Orszag- 
Tang vortex made up of a magnetic X-point centered on a stagnation point of the 
velocity. A recently developed spectral-element adaptive refinement incompressible 
magnctohydrodynamic (MHD) code is applied to simulate this problem. The MHD 
solver is explicit, and uses the Elsasser formulation on high-order elements. It 
automatically takes advantage of the adaptive grid mechanics that have been described 
elsewhere in the fluid context [Rosenberg, Fournier, Fischer, Pouquet, J. Comp. Phys. 
215, 59-80 (2006)]; the code allows both statically refined and dynamically refined 
grids. Tests of the algorithm using analytic solutions are described, and comparisons 
of the Orszag-Tang solutions with pseudo-spectral computations are performed. We 
demonstrate for moderate Reynolds numbers that the algorithms using both static and 
refined grids reproduce the pseudo-spectral solutions quite well. We show that low- 
order truncation-even with a comparable number of global degrees of freedom-fails 
to correctly model some strong (sup-norm) quantities in this problem, even though it 
satisfies adequately the weak (integrated) balance diagnostics. 
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1. Introduction 

In geophysical and astrophysical flows, the Reynolds numbers are large, and nonlinear 
terms appearing in the primitive equations lead to strong mode coupling and multiple 
scale interactions. Moreover, magnetic fields are often in quasi-equipartition with the 
velocity, as is the case in the Solar Wind or in the interstellar medium. However, direct 
numerical simulations (DNS) using regular grids cannot, even to this day, deal with 
such large Reynolds numbers R v . Doubling the grid resolution (and thus multiplying 
the Reynolds number by roughly a factor two) comes at a cost of increase in needed 
computer time by a factor of sixteen in three dimensions, assuming the temporal scheme 
is explicit; even when taking Moore's law into account, such an increase in R v can only 
be achieved roughly every six years. Thus, one is led to resort to more sophisticated 
techniques, such as, for example, turbulence modeling (see e.g. [21] ). However, in 
the case of coupling to a magnetic field, and using the magnetohydrodynamic (MHD) 
approximation valid for the description of the large-scale dynamics, few such techniques 
have been developed and tested (see however [26], [22l EHJ [15] ) - Another venue is to 
develop adaptive mesh refinement (AMR) methods. In this context, we examine in this 
paper the accuracy of an AMR code using spectral elements by comparing its output 
to exact solutions in simplified cases and to computations using a pseudo-spectral code 
at the same Reynolds numbers on a classical configuration of magnetic reconnection in 
two-dimensional geometry. 

We set up the equations in the next section [21 and in section [3j describe the 
numerical method for MHD developed within the context of the adaptive spectral 
element method presented in [30] for the two-dimensional Burgers equation. This section 
also presents test results for the method. In section H] we apply the method to the 
Orszag-Tang configuration, and compare with pseudo-spectral results. We consider 
effects of low order versus high order local approximations in section [51 and section [6] is 
the conclusion, in which we summarize the results, and offer some observations about 
the performance of the method and some directions for future work. 

2. Setup and theory 

2.1. Equations, code, and simulations 

For an incompressible fluid with constant mass density po, the magnetohydrodynamic 
(MHD) equations read: 



dru + u ■ Vu 



Vp + j x b + z/V 2 u, 



(1) 




(2) 



V-u = 0, V-6 = 



(3) 



where u and b are the velocity and magnetic field (in Alfven velocity units, b = B / ^/p^opo 
with B the induction and /iq the permeability); p is the pressure divided by the mass 
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density, and v and r) are the kinematic viscosity and the magnetic resistivity. The mode 
with the largest wavevector in the Fourier transform of u at initial time is ko- We also 
define the viscous dissipation wavenumber as k v = (e/V 3 ) 1 / 4 , where e ~ Uq/L is the 
energy dissipation rate, with Uq the r.m.s. velocity and Lq the integral length scale (see 
below). The Kolmogorov scale at which dissipation sets in is defined as Id = 2ff/k 1/ ; 
the expression for k v is based on a kinetic energy spectrum Ey(k) ~ e 2 / 3 fc~ 5 / 3 . A 
large separation between the two scales (k$ 1 ^> k^ 1 ) is required for the flow to reach a 
turbulent state with significant nonlinear interactions. 

In the absence of external forcing, viscosity and magnetic resistivity, the MHD 
equations in two space dimensions (2D) conserves the total energy: 

E= 1 -j(u 2 + b 2 )dx 2 , (4) 

together with the cross helicity H c = ~ J u -b <ix 2 , and the L 2 norm of the magnetic 
potential M a = | J a 2 dx 2 , with b = V x a. 

The Reynolds number is defined as R v = UqLq/v, where the integral length scale of 
the flow is defined as 

jEy{k)k^dk 

L ° = ^ jEy(k)dk ■ (5) 

The large scale turnover time can then be defined as t^l = Uo/L . We can also introduce 
the Taylor based Reynolds number R\ = UX/u, where the Taylor length scale A is given 
by 



jE v (k)k 2 dk / 

Length scales built with b and its energy spectrum Euik) can also be defined; the 
magnetic Reynolds number is R m = U L /r). 



3. Algorithm description for MHD 

For this work, we use a spectral element method [28], encapsulated within the 
Geophysical /Astrophysical Spectral-element Adaptive Refinement (GASpAR) code. 
Aspects of this code-in particular those regarding the dynamic grid refinement-have 
been described elsewhere [3D], where results were presented for the multi-dimensional 
Burgers (advection-diffusion) equation. Here, we extend that development in several 
distinct ways in order to solve (0Q)-(l3]), specifically, by adding the pressure term and 
the Lorenz force in the momentum equation and by taking into account the magnetic 
induction equation, and the divergence-free conditions on velocity and magnetic field. 
We solve equations ([I])-© in Elsasser form [6]: 

d t Z± + Z T ■ VZ ± + Vp- z^V 2 ^ - v T V 2 Z T = (7) 
V • Z± = , (8) 
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where 

Z ± = u±b 

and 

v ± = -{v±rf). 

Thus, we solve for u and b in terms of Z . 

Equations ©-([8]) can be recast into an equivalent variational form on domain D 
by defining the following function spaces: 

U 5 :=|w = ^^e M | G U X (D) V/i & w = w on dD^j (9) 

(10) 

H\D)-.= {f \feC 2 (D)kd x ,feC 2 (D)^}r, (11) 

where w = u, b. Let Z and p and their test functions, and g be restricted to 
finite-dimensional subspaces of these spaces: 

Z± G U N = U t fl P N (12) 
<f± G < = U f| P N (13) 
MeY M = D f| P^. (14) 

The basis for the velocity expansion in P^v is the set of Lagrange interpolating 
polynomials on the Gauss-Lobatto-Legendre (GL) quadrature nodes, and the basis for 
the pressure is the set of Lagrange interpolants on the Gauss-Legendre (G) quadrature 
nodes. Then, the equations (|7|)-((8|) can be written in weak form as [9] 

(^,d t Z ± ) GL + (( ± ,C T Z ± ) GL --(p, v-c^g = -^(vc^, vi±) GL (i5) 

Po 

(g,V-Z ± ) G = 0, (16) 

where C 1 * 1 := Z ± • V is the continuous advection operator, (•, -)gl, represents the inner 
product using quadrature on the GL nodes, and (•, -)g indicates inner product using 
quadrature on the G nodes. Thus, we use a staggered grid, where the quantities (u, b, 
Z ± ), on the GL nodes are continuous, while those on the G nodes (p) are not. This 
so-called P N — P N ^ 2 wa s chosen to prevent spurious pressure modes [191 E] • 

In the spectral element method, functions in and Y N ~ 2 are represented as 
expansions in terms of tensor products of basis functions within each subdomain, or 
element, the non-overlapping union of which comprises the domain [28J : D = Ufc=i E fe- 

By expanding Z^ and p in terms of their basis functions, inserting these expansions 
into (fT5l) - (|T6l) . and using the appropriate quadrature rules, we arrive at a set of semi- 
discrete equations in terms of spectral element operators: 
dZ ± 

M ~dt = ~ MCTZ t + Vjp* ~ u ± LZ t ~ u t LZ J (17) 
D j Zf =0, (18) 
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for the j th components of momentum. In this equation, M, L, and C, are the well 
known mass matrix, weak Laplacian and advection operators, respectively (e.g. , [30] 



and references therein) . The variables Z represent values of the Z^ collocated at the 
GL node points, while p^ 1 are values of the pressures at the G node points. Similarly, 
we denote the collocated test functions by q. The Stokes operators, Dj, arise from 
the quadrature rule in (I16p . in which the GL basis function and its derivative must be 
interpolated to the G node points, and multiplied by the G quadrature weights. In two 
dimensions, 

K 

(q, V • Z) G = £ (q k ) T (D k 1 Z k l + D k 2 Z k 2 ) (19) 

fc=i 

->±,k , 

for Z = Z . For regular rectangular elements, of lengths L*, 



where 
and 



d;=[^i®d, d*=(^)d®i, 



Uj = Vi^j&i 



n - n —1\ 

dr 4 



are, respectively, the weighted one-dimensional (ID) interpolation matrix mapping GL 
points to the G points, and the weighted ID GL derivative matrix interpolated to the 
G points [9] . In these definitions cr, are the weights corresponding to the G node points. 
Just like with the mass and Laplacian operators, the Stokes operators are written as 
tensor products of ID operators. 

The above equations, (1TT T) -( TT8T) . are correct for a single element, but they are not 
complete when we have multiple subdomains. In this case, we must ensure that all 
quantities in XJ N remain continuous across element interfaces. The manner in which 
this is done for advection-diffusion on (non) conforming elements was described in [30J. 
Using the Boolean scatter matrix, A c , the interpolation matrix from global to local 
degrees of freedom, and the masking matrix (that enforces homogeneous boundary 
conditions), \~\ , that were presented there, we find that we must replace the local Stokes 
operators above with 

Dj -> D Lj $A c n, 

where Df,j = diag fe (D ), and the D k are the matrices from above. In much of what 
follows, we will continue to use the local form of the Stokes operators, and simply 
recognize that the multiple subdomain form can be imposed. 

Note in (TTTT) the presence of different pressures for Z ± . As we show below, we will 
maintain the divergence constraints by solving for the pressures for both Z^. While 
analytically these pressures are the same, they serve as Lagrange multipliers [7] for 
their respective fields, Z ± . Given that each field has its own constraint that is solved 
independently of the other, numerically they will in general be different. 
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3.1. Time stepping 

Because we want to resolve all time scales (as well as spatial scales), we choose an explicit 
method for integrating (fT7|) - (fT%|) in time. The method we use is the m^-order Runge- 
Kutta (RK), and since the right-hand side of the equations has no explicit dependence 
on time, we can use the following formulation [U p. 109] in order to solve dU/dt = F(U): 

Set 

U = U n 
For k — m, 1, —1 

U = U n + \AtF(U) 

End for 

U n+l = U. 

Considering one component of the Elsasser variables and following the above RK 
algorithm, we write for each iteration (recall (fT7j) ) 

Zf = Zf n - M~ 1 (MC T Z J ± - D]p ± + v±LZf + u T LZj). (20) 

We insist that each RK stage obey (1181) in its discrete form, so multiplying ( l20il by 
Dj, summing, and setting the term D J -Z^ = leads to the following pseudo-Poisson 
equation for the pressures, p ± : 

D^M^Djp* = D j gf, (21) 

where the remaining inhomogeneity 

gf = D.M^MC^ + u ± LZf + u T LZj) - DjZf. 

The pseudo-Laplacian operator on the left-hand-side of (l2~Tj) also arises from a second 
order implicit time descretization of the spatially discretized equations ffTTj) - ffT8l) as 
shown in P,[8]. In this formulation, even if Z^ is not divergence-free, the partial update 
after this RK stage will be. The inverse mass operator, M _1 , must be computed for a 
given grid configuration. For conforming elements, this matrix is trivially inverted since 
it is diagonal. For nonconforming discretizations, the mass matrix is lumped in order 
to recover a diagonal matrix that can be inverted straightforwardly [8]. Equation (|2~T1) 
is solved using a preconditioned conjugate gradient method (PCG; see [33j [36] , and also 
[30] for modifications required for nonconforming elements). For preconditioning, we 
use a block Jacobi preconditioner computed using either a fast diagonalization method 
[5], or a direct inversion. 

Thus, at each timestep, m RK stages are computed, and each stage solves fl2~Tl) 
twice, once for Z + , and once for Z~ leading to 2m pseudo-Poisson solves at each time 
step. The MHD solver is encapsulated within a class as are all solvers in the code. 

Currently the solver considers the Elsasser variables, Z , to be auxiliary variables, 
so a transformation is done from the native u, b, and back. However, it is easy to 
create a new constructor for the solver object, so that Z^ 1 are the primary quantities of 
interest. In all results presented in this work, we choose m = 2 for the RK order. 
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Figure 1. (left) Elemental decomposition and streamlines of the steady-state 
Kovasznay solution as determined by time marching for a case where p — 6. (right) £2 
error vs. polynomial degree, showing spectral convergence of the Kovasznay solution. 



3.2. Tests with exact solutions 

To test the algorithm described above, we first consider two test problems that have 
analytic solutions. The first case tests the algorithm without magnetic field (b = 0) 
field-ie., solves the Navier-Stokes equations-and the second case tests the full MHD 
algorithm. 

For Navier-Stokes, we use a steady state solution of a (Kovasznay) flow field behind 
a periodic array of cylinders p2] (see also [5]): 

u x = 1 — e Xx cos 2iry 

u y = -^e Xx sin 2iry, 
Ztt 

where A = -y — \J ^ + 4n 2 . We initialize the grid with the solution, and march to 
a steady state (typically to t m 2), where we compare the numerical solution with 
the analytic solution. For the following test, the grid we use is conforming, and non- 
adaptive (see figured]); the time step is fixed at 1 x 10~ 3 , while R v = \ jv — 40. Dirichlet 
boundary conditions are set from the analytic solution. We plot the £2 norm of the error 
as a function of polynomial degree, p in figured] (right), in order to show convergence. 
Clearly, the solution is converging spectrally as desired. 

The next test looks at a steady Hartmann flow, consisting of a flow of viscous 
conducting fluid between two parallel plates, with separation 2a, and with a magnetic 
field, B , applied in a direction perpendicular to the plates. The solution is [T8] : 
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Figure 2. £2 error vs. polynomial degree p, showing spectral convergence of the 
Hartmann solution for u x and b x . Similar convergence is seen for u y . 

y/a smhH a — smh(yH a /a) 
cosh H a — 1 
where the Hartmann number, 

a aB 
^Jvr\ 

and 5 represents a boundary layer thickness over which the velocity compared with that 
in the central region decreases at the top and bottom plates. 

Again, we initialize a grid of K = 4 x 2-elements-whose lower left and upper 
right global grid vertices are (0, 0), and (4, 2)-with the analytic solution, at fixed p and 
march to steady-state, comparing the numerical an analytic solutions. For these runs 
we use a Courant-limited timestep. Dirichlet boundary conditions are again set from 
the analytic solution. To drive the flow, we add a constant x- momentum forcing, ~{ x , to 
the right-hand-side of §\§ (or, equivalently, to ©) such that 

ul H a smhH a 
a R v cosh H a — 1 

For all our tests, we set H a = 4, R v = ^ = 40, u = 1, and B — 1; the choice of 
H a > 1 suggests a flow where the magnetic field is dynamically significant. In figure [2] 
we present our convergence results. As expected, we again see spectral convergence. 

4. The Orszag-Tang vortex 

The Orszag-Tang vortex (hereafter, OT; [27]) is a simple configuration with a magnetic 
X-point centered at a stagnation point of the velocity. It is best expressed in terms of 
the stream function \I/ (with u = V x ^ z) and the magnetic potential. It reads 

$ = 2a [cos(27rx) + cos(27n/)] 

a z = 2cos(27rx) + cos(2ttc/). 
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The total energy Et = Ey + Em, for a = 1, is equally divided initially between its 
kinetic and magnetic components Ey and Em, both equal to 2, and the initial correlation 
coefficient p c is equal to 41%, where p c = J^ 2 . This configuration is known to develop 
several current sheets in a time of order unity in these units; such current sheets further 
destabilize in time through a tearing mode instability embedded within a turbulent flow 
and develop numerous complex small-scale structures (see e.g. [29]). 

In this section, we apply the algorithm described in Section 3 to the OT problem. 
For adaptive runs, a spectral estimator refinement criterion is used [161 1201 EH] that 
estimates the solution error. If, in a given element, this error is greater than a specified 
tolerance, e est , the element is tagged for refinement. The nominal resolution of the 
adaptive runs is given by an equivalent resolution, which is computed by 

N cq = pN 2 £ ^, 

where N Q is the initial number of elements in either direction (the base grid), and £ max 
is the maximum refinement level. All computations are performed on a periodic grid of 
dimension [0, l] 2 . We compare the spectral element solutions with those obtained from a 
well-characterized pseudo-spectral code that has been used to produce numerous results 
cited in the literature [131 El • 

The total dissipation in the flow is defined as 

V T = -v < u 2 > -7] < f >, (22) 

where u = V x u is the vorticity and j = V x b is the current density. It is a global 
quantity, as is the total energy Et, and is characteristic of the dynamical evolution 
of the flow as a whole: as small-scale gradients in both the velocity and the magnetic 
field develop through non-linear interactions, current and vorticity sheets form and 
dissipation sets in at a time of order unity. The temporal evolution of Et and V>t are 
shown in figure [3] for the GASpAR MHD run (solid line) and the pseudo-spectral code 
(dotted line) for a fiducial run with R v = 10 3 7r; the initial grid is K = 8 x 8 with £ max = 3 
and p = 8 (implying N eq = 512), and e cst = 1 x 1CT 5 . The adaptive code reproduces the 
temporal characteristic times (as well as secondary maxima in T>t at lower Reynolds 
number, shown in figure [TUl in the context of a fixed grid, see below); it also reproduces 
the amplitude of the global phenomenon of reconnection of magnetic field lines and 
ensuing dissipation of energy. 

The total number of degrees of freedom (d.o.f.) normalized by the number of modes 
in the pseudo-spectral run used in this comparison during the temporal evolution of the 
flow is also shown in figure [3] in black. The d.o.f. quickly increases because the initial 
grid was coarse (K = 8x8 elements, at p = 8). After this, the d.o.f. grows regularly 
in anticipation of the peak in total dissipation at a slightly later time; note that it is 
about one-half that of the pseudo-spectral computation. Other runs at lower Reynolds 
numbers indicate that, beyond this peak, the number of d.o.f. is roughly constant until 
about t = 1.8, where the grid appears to anticipate a secondary peak in T>t that begins 
at about t = 2.5, and corresponds to renewed reconnection of current layers [29] (not 
shown) . 
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Figure 3. (a) Comparison of enstrophy (red curve) with a pseudo-spectral code with 
512 2 d.o.f. (black curve), as a function of time. The curves overlie each other, (b) 
Degrees of freedom in the adaptive run normalized by the total number of d.o.f. in the 
pseudo-spectral run, as a function of time, for the refinement criteria based on u x and 
B x (black), and j and u> (green). 

As the relative number of d.o.f. increases toward unity, the AMR becomes less 
efficient. In this run, the grid refines on the native solver variables u and b, and by 
the time the nonlinear regime begins, the small scales dominate the flow requiring finer 
elements to remain on the grid to resolve them. The number of d.o.f. also depends 
on the criterion of refinement and on e cst ; note that here e es t is set tightly so that the 
grid is more likely to refine. It will be left to future work to investigate this and other 
refinement criteria more systematically, in order to determine what range of parameters, 
and, indeed, which criteria, are the most robust and efficient for adaptive modeling of 
flows such as OT. We have freedom in the code also to define new variables to which to 
apply our refinement criteria, providing yet another avenue for investigation (see also 
the discussion below around figure [6]). As an illustration, we also show in figure [3] the 
number of d.o.f. for a spectral criterion of refinement, based this time on the vorticity 
and the current. In this case, the refinement of the grid starts at a later time, once 
the strong gradients have developed, and thus the run is roughly twice as fast; however, 
once we approach the saturation of the growth of small scale production, the number of 
d.o.f. for the two refinement criteria are comparable. 

Conservation of energy (and that of the other quadratic invariants) is an essential 
feature in the detailed dynamical evolution of a turbulent flow; it is at the foundation 
of the concept of energy (or invariant) cascade and leads, assuming a constant energy 
dissipation rate e within the cascade, to the celebrated Kolmogorov energy spectrum 
E v (k) ~ k~ 5/3 . Note that in MHD, the power law followed by the energy spectrum 
in the inertial range is less clear, and other spectra can be postulated a priori on the 
basis of Alfven wave propagation, including an anisotropic component of the spectrum 
linked with the bi-dimensionalization of the flow in the presence of a strong uniform 



Spectral adaptive mesh refinement in MED 



11 




k 



Figure 4. Energy spectra for the fiducial run at a time near the peak in enstrophy. 
The black curves represent the spectra computed from the pseudo-spectral method, 
and the red, from GASpAR. The upper curves are the magnetic energy, while the lower 
are the kinetic energy. 

magnetic field. Such power laws are barely observable, due to a variety of reasons. In 
the laboratory or in observations of geophysical flows, the instrumentation has cut-off 
frequencies, and in numerical simulations, the resolution is barely sufficient to be able 
to resolve the ranges necessary to follow the evolution of the flow, namely the energy- 
containing range around k ~ ko, the inertial range where loss- less transfer of energy 
occurs, and the dissipation range. For example, for fluids in three dimensions, it was 
shown in [l] that the Kolmogorov energy spectrum, on resolutions of up to 1024 3 points, 
occurs on a small range of wavenumbers (slightly more than two octaves) and is followed 
by a shallower power law, named a bottleneck effect. 

The energy spectra close to the maximum of enstrophy, at t m 0.85, are given 
in figure HI the spectra for the AMR run are computed on an irregular grid using the 
algorithm derived in [10]. We see that the agreement is quite good between the adaptive 
spectral element run and the pseudo-spectral run. It is interesting to note that, while the 
pseudo-spectral case uses dealiasing, no explicit dealiasing is required in the GASpAR 
run. 

A further test of the code is to verify that the nonlinear terms of the primitive 
equations do preserve the invariants, as spectral methods are expected to be 
conservative. This can be done by computing the time derivative of the invariant (say, 
using an algorithm of order 2), and comparing it to its theoretical value. In the case of 
the total energy, for example, the dissipation T>t given in equation (I2~2l is the theoretical 
value. In figure [5] we show the degree to which the energy is preserved and observe again 
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time 

Figure 5. Time variation of the energy invariant compared to its exact dissipation 
when v ^ 0, J] ^ 0. The blue curve is the time derivative and the red, the dissipation 
term. Both curves are normalized by the maximum in the dissipation. There is no 
discernible difference between the two curves. 
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Figure 6. Contours of current density, j, for the fiducial run, showing the refined 
grid. Each yellow box is an element; the grid lines indicating the node points within 
each element are not shown. 

an excellent agreement. 

A snapshot of the current density shown in figure [6] near the peak in enstrophy 
indicates that, as expected, the grid refines in the region of the current sheets (and 
vorticity sheets, which are known to be almost co-located). However, when compared 
with runs from other authors [H] using finite difference methods and configuration space 
refinement on the gradients of the basic variables, there appears in our run to be more 




refinement than is initially expected outside of the sheets. In the case of the spectral 
refinement criterion applied to u and b, which are native to the solver, the criterion 
appears to capture the variation in the curvature of u and b near the current sheets to 
an extent that may not be required since, with the same e es t, half the d.o.f. is needed 
at early times with the spectral criterion based on j and u. 

We now examine the temporal evolution of the current maximum in this fiducial 
run, as shown in figure H As is well known, the development of singularities in such 
flows may be diagnosed by a l/[£* — t] behavior [3], with t* the time of singularity. 
As such, the temporal evolution of extrema is an important feature of turbulent flows; 
from a physical point of view, such extreme events are also the locus of anomalous 
diffusion and concentration of (say, chemical) tracers so that a faithful reproduction of 
such events may be required. The linear (exponential) growth of the current instability 
is well reproduced but as we enter the nonlinear phase, studied in the literature in the 
context of the nonlinear growth of the tearing mode instability [32] , errors appear when 
using the AMR codes, errors that are comparable for the two refinement criteria just 
described (see figure CJ) . As already shown in three dimensions [21] , the nonlinear phase 
is highly nonlocal, which indicates that many scales are interacting in the development 
of current instabilities. This discrepancy (not observable in the £2 norm) might be 
remedied by tightening the refinement criteria. However, the accuracy of the code, 
as measured primarily by the polynomial order in each element for a fixed number of 
elements and without refining the grid, is a parameter that must also be examined, and 
this is done in the next section. 
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Table 1. Parameters used in the simulations described in Section 5; p is the polynomial 
order in each direction for each element, N eq is the linear grid resolution, and E is the 
number of elements in each coordinate direction such that the total number of elements 
is K = Ex E. 
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5. High versus low order 

We now consider the behavior of the OT solutions when the number of global degrees of 
freedom is kept (roughly) constant, while the polynomial truncation (degree) varies. In 
the first series of runs, we have v = rj = 0.025 (R v = 807r), and these runs are compared 
with a pseudo-spectral run with resolution of 128 2 grid points. Here, as stated before, 
we use a static conforming mesh for each run so that the results will not be affected by 
refinement and coarsening criteria, since we want to focus on the effect the order of the 
method has on the results. In table [U we present the relevant run parameters. 

In figure El is shown the profile of J max as a function of time for each of the iV eq ^128 
runs. It is seen immediately that the low order truncation does not yield accurate values 
for J max , particularly as the peak in total dissipation is reached near t = 1.1. Clearly, 
Jmax converges to the correct solution as p increases (compare the red curve for p = 5 
and the black curve for the pseudo-spectral run). The maximum of vorticity behaves in 
much the same way as J max . 

At this stage, it is desirable to compare the accuracy for a pseudo-spectral code 
with N points per linear dimension, e ps , to the error, e se , of a spectral element code with 
E elements per linear dimension, each with polynomials of order p. Omitting prefactors 
with slower variations in the truncation orders, we have for the pseudo-spectral error 
[H p. 400] that 

e ps ~ Ax N ~ 1/N N 
and for the spectral element error bound (5j p. 273], that 

e se ~ h mill{p ' s) p~ s 

where h ~ 1/E is the uniform element length, and s is the smoothness of the exact 
solution. It is clear that in practice s < p, since the derivatives for s > p cannot be 
computed. We thus choose s = p so that the function is sufficiently smooth to allow 
spectral convergence. 

Equating the logarithmic errors immediately shows the relationship between N, E 
and p, namely: 

iVlogiV~ plog(p E). 

The above scaling argument allows for choosing a range of values of polynomial 
order under simple assumptions. Let us say the Reynolds number is doubled from a 
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Figure 8. Maximum of current J max as a function of time for the first series of 
N cq w 128 runs in Table [TJ The black curve is the pseudo-spectral case; red is the 
p = 5 run; dark blue is the p = 4 run; green is the p = 3 case. Note that the p = 5 and 
pseudo-spectral results are nearly coincident. 

well-adjusted run; roughly speaking, we need to double the resolution of the pseudo- 
spectral code from N to 2N grid points per linear dimension. Let us also assume 
that we double the number of elements in the spectral element code. Then, under the 
reasonable assumption of large N, we find that the polynomial order needs to be doubled 
as well. Though empirical, this criterion does indicate that the polynomial order needs 
to increase with Reynolds number. It should be noted that, whereas the pseudo-spectral 
code, in the preceding example of figure El uses a truncation of order 128, the equivalent 
spectral element code uses polynomials of order 5, substantially smaller, in order to 
achieve comparable accuracy. 

As a specific example, an examination of figure [8] indicates that, with E = 26, 
setting p = 5 leads to a satisfactory computation of J max - Let us now double the 
resolution of the pseudo-spectral code to N' = 256 and take E' = 52. As can be seen in 
figure [9J it is indeed the case that p > 8 gives an accurate representation of the dynamics 
of the flow at that enhanced Reynolds number. According to the scaling relationship, 
if we were to retain p = 5, we would need E ~ 10000 in order to reach the same level 
of accuracy at the higher R v . For completeness, we show in figure [10] the corresponding 
£.2 norms for the total energy (top) and the total generalized enstrophy < u) 2 + j 2 > for 
the N eq ps 256 cases; the different runs cannot be discerned. 

6. Discussion and Conclusion 

We have presented an explicit spectral element method for solving the equations of 
incompressible magnetohydrodynamics. The method is developed within the context 
of an existing spectral element code (GASpAR) that provides many of the spectral 
element operators required of the MHD algorithm, and that also offers an adaptive, 
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Figure 9. Maximum of current J m ax as a function of time for the second series of 
N eq ~ 256 runs in Table [TJ The black curve is the pseudo-spectral case; red is the 
p = 8 run; dark blue is the p = 6 run; green is the p = 4 case; cyan is the p = 3 run. 
Note that this time the p = 8 and pseudo-spectral results are nearly coincident. 
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Figure 10. Plots of total energy (top) and enstrophy (bottom) for the N eq « 256 runs 
in Table [TJ The color scheme is the same as in figure [5J but the curves lie on top of 
one another. 

nonconforming mesh algorithm. The new operators that arise in the explicit MHD 
treatment, have been defined. Here, we have described tests that compare the 
numerical results with analytic solutions and establish that the method achieves spectral 
convergence in the case of conforming elements (for preliminary tests on a static 
reconnection problem, see [31]). We have then applied the method to a challenging 
problem in the literature, the so-called Orszag-Tang flow, which allows for magnetic 
reconnection and the development of current sheets. The OT runs are compared with 
well-tested pseudo-spectral solutions as a baseline, and found to agree well. We find 



Spectral adaptive mesh refinement in MHD 



17 



that the quadratic £2 diagnostic quantities are insensitive to variations in polynomial 
degree, while the sup-norm quantities, such as the maximum current density, are not 
accurate at low order truncations. Because such sup-norm quantities are the foundation 
for criteria of the development (or not) of singularities in Navier-Stokes and MHD flows 
[21 [3], it may be of some importance to be able to solve for them accurately. 

It is worth pointing out that the spectral element method described in this paper 
is more costly in terms of computational time than the pseudo-spectral method with 
which we compare our results, the latter being optimal for periodic boundary conditions. 
This is despite the fact that the spectral element method requires only nearest-neighbor 
communication, while the pseudo-spectral method requires all-to-all exchange of data 
in the fast Fourier transform algorithm. This performance issue can be traced directly 
to the solutions of the pseudo-Poisson equation ( T2TT) that are required in order to 
maintain the divergence constraints. The pseudo-Laplacian operator is known to be 
ill-conditioned, primarily because of a one-dimensional null space. Indeed, for the 
iVgq = 512 2 equivalent runs, we see typical iteration counts for each RK stage of 
~ 250, which further increase-albeit reasonably slowly-in the OT runs as the enstrophy 
maximum is approached; we see no significant reduction in the PCG iteration count 
when we attempt to remove this null-space, but more work is needed in this area. [It 
should be noted that the scaling of the spectral element code on multiprocessors is very 
good, but we refer here to single processor performance.] Furthermore, when dealing 
with incompressible MHD fluids at low magnetic Prandtl number, as encountered in the 
laboratory and in the liquid core of the Earth, there is a need for accurate simulations 
of the generation of magnetic fields in turbulent flows with complex boundaries, in 
conjunction with several ongoing experiments (see e.g., (35J, [121 [331 [25] ) • Spectral element 
codes, which encompass easily a variety of boundary conditions and geometries may be 
useful in this context. 

There are a number of ways to speed up the spectral element solutions. One is 
to implement a more sophisticated preconditioner, and we are making progress in this 
regard that will be reported on elsewhere. Another is to relax the degree to which the 
divergence constraints are maintained, by increasing the PCG convergence tolerance to a 
less stringent value. Still another, as alluded to in section HI is to find refinement criteria 
that optimize the number of degrees of freedom for a desired quadrature or truncation 
error. We have seen that the adaptive mesh code can provide a substantial savings 
in work because the number of degrees of freedom, as shown in several instances (e.g. 
[30]). can be reduced. The code is presently being extended to include compressibility in 
view of the many MHD applications we have in mind in the astrophysical context (solar 
wind, magnetosphere, solar convection zone and corona), in which case the performance 
problems associated with this divergence-free constraint should be alleviated. 
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